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Hard constraints imposed in statistical mechanics models can lead to interesting thermodynamical 
behaviors, but may at the same time raise obstructions in the thoroughfare to thermal equilibration. 
Here we study a variant of Baxter's 3-color model in which local interactions and defects are included, 
and discuss its connection to triangular arrays of Josephson junctions of superconductors with broken 
time-reversal symmetry and kagome networks of superconducting wires. The model is equivalent 
to an Ising model in a hexagonal lattice with the additional constraint that the magnetization of 
each hexagon is ±6 or 0. Defects in the superconducting models correspond to violations of this 
constraint, and include fractional and integer vortices, as well as open strings within 2-color loops. 
In the absence of defects, and for ferromagnetic interactions, we find that the system is critical for 
a range of temperatures (critical line) that terminates when it undergoes an exotic first order phase 
transition with a jump from a zero magnetization state into the fully magnetized state at finite 
temperature. Dynamically, however, we find that the system becomes frozen into domains. The 
domain walls are made of perfectly straight segments, and domain growth appears frozen within the 
time scales studied with Monte Carlo simulations, with the system trapped into a "polycrystalline" 
phase. This dynamical obstruction has its origin in the topology of the allowed reconfigurations in 
phase space, which consist of updates of closed loops of spins. Only an extreme rare event dominated 
proliferation of confined defects may overcome this obstruction, at much longer time scales. Also as 
a consequence of the dynamical obstruction, there exists a dynamical temperature, lower than the 
(avoided) static critical temperature, at which the system is seen to jump from a "supercooled liquid" 
to the "polycrystalline" phase within our Monte Carlo time scale. In contrast, for antiferromagnetic 
interactions, we argue that the system orders for infinitesimal coupling because of the constraint, 
and we observe no interesting dynamical effects. 

PACS numbers: 00000 



I. INTRODUCTION 

Systems with hard constraints often display interesting 
thermodynamic properties such as infinite order phase 
transitions or, on the contrary, very sharp first order 
phase transitions. Many of these models can be de- 
scribed in terms of vertex models and some of them are 
exactly solvable. Examples of such systems are given by 
dimer models 1 , the planar ice model 2 or the three color- 
ing model of the hexagonal lattice^. 

It is very natural to ask whether the hard constraint, 
which leads to the interesting thermodynamics, may at 
the same time pose obstructions in the (possible?) path 
to thermal equilibration. In essence, equilibrium proper- 
ties require averages over all the configurations allowed 
by the constraint, weighted in accordance with the ap- 
propriate Boltzmann-Gibbs distribution. Dynamically, 
the system must sample the different allowed states in 
a manner that satisfies detailed balance. However, leap- 
ing from an allowed configuration to another might re- 
quire large rearrangements, and physically one must in- 
vestigate which mechanisms could possibly lead to these 
moves in phase space and what are the corresponding 
time scales. Sometimes the constraint forbids any lo- 
cal rearrangement of the system (as in the present case), 
and it ought to be softened in order to allow for a lo- 
cal dynamics. The system then evolves by formation of 



constraint-violating defects that propagate and recom- 
bine. 

Plenty of issues arise regarding the dynamical genera- 
tion and recombination of defects, which depend on the 
microscopic details of the physical system, and the en- 
ergetics of the states outside the manifold of constraint- 
satisfying states. For example, paying the energy cost to 
create a defect already slows down the dynamics; how- 
ever, this waiting for the defect generation simply rescales 
the time scales for dynamical evolution in a trivial way. 
More interesting are those issues related to the possible 
energy costs for moving defects around. In particular, 
if the microscopies are such that the defects (when cre- 
ated in pairs) are confined, one would expect further and 
non-trivial slowing down of the dynamics. 

Glassy behavior in constrained 3-color models with in- 
finite range interactions has indeed been recently found 
by Chakraborty, Das, and Kondev 4 . This is an interest- 
ing example of glassy behavior in a Hamiltonian model 
without quenched disorder, where it was found that the 
characteristic time scales obeyed a Vogel-Fulcher law as 
the temperature approached a dynamical transition tem- 
perature, mimicking fragile structural glasses. In order 
to maneuver within the phase space of allowed states, 
non-local loop dynamics was implemented. 

In this paper, we study variations of the Baxter 3-color 
model with short range interactions and discuss the possi- 
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ble mechanism for defect motion. In particular, we argue 
that the loop updates used by Chakraborty et alX cor- 
respond to the unbinding of certain defect pairs that are 
deconfined, and thus they are the least costly mechanism 
for dynamical evolution. We find that finite range fer- 
romagnetic interactions lead to a frozen "polycrystal" , 
as opposed to a fragile glass as in the case of infinite 
range interactions. We present two possible experimen- 
tal realizations using lattice arrays of superconducting 
devices that could in principle be experimental settings 
for studying sluggish relaxation or non-equilibrium ef- 
fects in Hamiltonian systems without quenched disorder. 

In section |H] we present in detail the 3-color model, 
and show that it is equivalent to an Ising model on a 
hexagonal lattice, with the constraint that the magneti- 
zation of each hexagon must be ±6 or 0. In the Ising 
language the extra interaction that we add to the 3-color 
model has a simple form: it is a nearest neighbor spin- 
spin interaction. Such interaction is present in the pos- 
sible experimental realizations of the model in two dif- 
ferent 2-D superconducting geometries. Because of the 
constraint imposed on the plaquettes, the system is crit- 
ical in the absence of two-spin interactions (J = 0) and 
is described by a c — 2 Conformal Field Theory (CFT) 5 . 
In Sec. IIIII we use this description to argue about the 
behavior of the model in the presence of non-zero two- 
spin interactions. While for arbitrarily small antiferro- 
magnetic coupling (J < 0) the system orders, it remains 
critical for small ferromagnetic coupling (J > 0). The 
CFT description near the J = point is ill-suited for 
strong couplings. In this regime we use instead a Clus- 
ter Mean Field Method (CMFM) which has proven to be 
very accurate in describing constrained system such as 
the ice models. We find a strong first order phase tran- 
sition where the system jumps from the disordered con- 
figuration to the Fully Magnetized Ferromagnetic State 
(FMFS). 

When the hard constraint is softened, defects are al- 
lowed in the system at a high energy scale U, which enters 
in the defect formation energy and in the defect pair in- 
teractions. In Sec. lIVI we discuss the role of these defects 
and their implications in the dynamics of the system. 
In the superconducting realizations there are a number 
of different defects: fractional vortices, integer vortices, 
and open segments of closed two-color loops. Integer and 
fractional vortices can be shown to be confined below a 
Kosterlitz-Thouless transition temperature that can be 
rather high depending on the energy scale U. Thus, 
these defects are rather ineffective as a mechanism to 
move from one allowed state to another. We show, on 
the other hand, that the end points of open segments of 
closed loops made of two alternating colors are decon- 
fined, they can move around and travel a whole closed 
loop, and therefore they are the main actors for the evo- 
lution of the system. For defect formation rates much 
smaller than the defect recombination rates, this evolu- 
tion corresponds essentially to the loop dynamics that we 
use in the present paper. 



In Sec. [V] we study the dynamics of the constrained 
system. By fitting the value of the free energy for the 
disordered state as a function of temperature and com- 
paring it to the one of the ordered state we first obtain an 
accurate estimate for the transition temperature, which 
is in good agreement with the result from the CMFM. 
We then show that there is no sign of the above men- 
tioned thermodynamic transition to the FMFS. The sys- 
tem instead becomes supercooled and undergoes a lower- 
temperature non-equilibrium transition from the super- 
cooled liquid phase to a frozen "polycrystalline" phase. 
The transition shows features that are characteristic of 
first order phase transitions, such as a hysteretic behav- 
ior as a function of temperature. The underlying physics 
behind this phenomenon is understood by studying the 
spin-spin autocorrelation function as well as the evolu- 
tion of the internal energy and other physical quantities 
when we cool the system at different cooling rates or after 
a quench from infinite temperature. 



II. THE MODEL AND ITS POSSIBLE 
EXPERIMENTAL REALIZATIONS 

In this section we review Baxter's 3-color model, and 
present two of its possible experimental realizations in 
lattices of superconducting devices in some detail. We 
show that the 3-color model and these two realizations 
can be described as an Ising model on a hexagonal lat- 
tice, with a plaquette constraint of ±6,0 for the sum of 
the spins around each hexagon. It is important to no- 
tice that while the 3-color model is only Z2 symmetric 
in the Ising spin representation, the superconducting re- 
alizations have a larger Z2 x U(l) symmetry due to the 
superconducting phase. This difference is particularly 
relevant for the possible defects that can originate in an 
allowed configuration and for their dynamic behavior. 

The one extra ingredient that we add to Baxter's 3- 
color model is a local interaction. In the Ising spin rep- 
resentation, this interaction takes the form of a nearest 
neighbor spin-spin interaction. It has the effect, in the 
3-color model, of aligning or not bonds of the same color 
on neighboring sites. The extra interaction is responsible 
for all the interesting thermodynamical and dynamical 
effects that are studied in this paper. Moreover, in the 
lattices of superconducting devices these interactions are 
always present. 



A. The three-color model 

The three-color model consists of vertices having three 
bonds of different colors: A, B and C. These different 
colors can be thought of as three different phases differing 
pairwise by ±27r/3, which is how we will later connect 
the model to arrays of superconducting devices. One 
can naturally associate to each vertex a chirality spin ±1 
depending on the counterclockwise or clockwise ordering 
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of the phases, as shown in Fig. ^ A hexagonal lattice is 
constructed with these vertices by connecting the bonds, 
where the connected bonds must share the same color. 
As we show below, the chirality spins cannot adopt an 
arbitrary configuration. Indeed, the spins must satisfy 
the constraint that their sum around any hexagon of the 
lattice is ±6, 0. On the other hand, given an allowed 
configuration of the spins, there are clearly three different 
corresponding color configurations, since any global even 
permutation of the colors in the lattice gives rise to the 
same spin configuration. In the absence of any kind of 



tex belongs to one and only one of such loops. Thus, by 
simply removing all the bonds of one of the three colors 
(say C), we realize one of the three possible simultaneous 
mappings of the system to a fully packed loop configura- 
tion on the hexagonal lattice which, at large scales, can 
be described by an su(3) level 1 Wess-Zumino-Novikov- 
Witten (WZNW) mode£. 

The 3-color model becomes even richer when we intro- 
duce a nearest neighbor spin-spin interaction in the Ising 
representation, which we do in subsection III El after we 
discuss the experimental realizations right below. 




FIG. 1: The gluing of the ABC vertices gives Baxter's three 
coloring model on the hexagonal lattice. To every vertex we 
can associate a chirality spin depending on the order in which 
the three colors appear counterclockwise around the vertex: 
+/— for even/odd permutations of the sequence ABC. 

interaction this model corresponds to the Baxter's three 
coloring model on the hexagonal lattice. The partition 
function Z has a purely entropic origin and its value is 
given by the number of ways of coloring the bonds of 
the hexagonal lattice. This number is known to grow 
exponentially with the system size. Indeed, Baxter solved 
exactly this model and showed that Z = W N for large 
values of the number of sites N, where W = 1.2087... is 
the entropy per sita^. 

It is worth discussing in detail how the system can re- 
arrange from one allowed configuration to another. No 
single-bond flip or double-bond exchange is allowed with- 
out violating the constraint in the neighboring vertices. 
However, we can notice that by choosing one vertex and 
two colors, say A and B, we can uniquely define a loop 
by taking the sequence of ABAB... bonds starting from 
the chosen vertex. The loop must be non self-intersecting 
and closed, the last property holding only if the system 
has periodic boundary conditions. Clearly, if we pick one 
such loop and we flip the color sequence, say ABAB... to 
BAB A..., the color constraint is preserved. These loop 
flips (or updates) provide a mechanism for the system to 
move around the phase space of allowed configurations. 
In Sec. II VI we will show how the loop updates originate 
from local constraint-violating defects. 

Notice that, given any allowed configuration, every ver- 



B. The Josephson junction array of 
superconductors 

A possible experimental realization of the model is 
given by a Josephson junction array of triangles of a su- 
perconductor with broken time reversal symmetry. For 
example, there is experimental evidence of a p x ± ip y or- 
der parameter in the compound S^RuO^^; here the two 
possible states p x ± ip y correspond to the chirality spin 
±1 defined above. The same geometry we propose here 
with p±ip states has also been studied by Moore and Lee, 
who in addition to the p-wave states have also looked at 
d ± id superconductors^, believed to be realized by the 
recently discovered hydrated cobalt oxide compounds. In 
their work, they have also discussed other type of arrays 
in triangular and square lattices. 




FIG. 2: An example of the correspondence between the 
Josephson junction array and the three-color model, provided 
we identify the three colors with the values of the phases of 
the order parameter in the middle of each triangle edge. No- 
tice that ferromagnetic order among nearest neighboring spins 
corresponds to aligning the bonds with the same color along 
the same direction. 

In the p x ±ip y Josephson junction arrays, the three col- 
ors correspond to the three relative phases of the order 
parameter in the middle of each of the edges of the trian- 
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gles, which differ by ±2w/3 (see Fig. 01. (To be precise, 
the phase of the order parameters are defined in momen- 
tum space; but, as it can be deduced from the analysis 
carried out in Appendix [S] one can think in real space 
by considering the phases for the momenta that point 
along the directions perpendicular to the three faces of 
each triangle.) The superconducting order parameter of 
each triangle has also an overall U(\) degree of freedom. 
Therefore, at the center of each of its three edges, one 
can define a phase Q, La = 6^ ± ^-a for the triangle at 
site i, along its a-th edge (a = 0, 1,2), where the edges 
are labeled counterclockwise starting from the horizontal 
one (see Fig.[3J). The ± sign corresponds to the chirality 




FIG. 3: Labeling of the edges of the up and down triangles, 
with the relative unit vectors e; jCl . While the chirality spins 
a p sit at the centers of the triangles, the "gauge" fields (an 
example of which is shown in one of the triangles) sit at the 
midpoints of the segments joining the centers of the triangles 
to the corresponding edges. Examples of the (7(1) phase 8 r 
and of the edge phases Q q , a , a = 0, 1, 2 are also shown. 

<7j = ±1 of the p x ± ip y state at site i. The Josephson 
coupling —Ucos{9i ta — 9j_ a ) along an edge shared by two 
neighboring triangles tends to align the phases 6i_ a and 
8j_ a . In the U — > oo limit one recovers Baxter's 3-coloring 
model, modulo a global U(l) phase. Notice that, in this 
infinite U coupling limit, the only difference between this 
system and the 3-color model (in the spin representation) 
described in the previous subsection is a Z 2 x U(l) sym- 
metry instead of a simple Z 2 symmetry. We will show 
in Sec. IIVI how this difference allows for a wider variety 
of defects in the Josephson junction array rather than in 
the three coloring model. 



C. The kagome network of superconducting wires 

Another (related) realization of the 3-color model is 
given by a superconducting kagome wire network in the 
presence of a magnetic fieldi2*ii*i£ such that the magnetic 
flux per triangular plaquette is one half of a flux quantum 
(/ = 1/2). Using a Ginzburg-Landau analysis, Park and 



Husei^ showed that the possible superconducting phases 
must have a gauge-invariant phase change around each 
elementary triangle equal to ±7r, and a gauge invariant 
phase change along each wire segment equal to ±7r/3. 
They also show that the allowed minimum free energy 
states of this model are equivalent to ground states of 
the XY kagome antiferromagnet, which are in one-to-one 
correspondence to the three-color model configurations, 
modulo aC/(l) phase analogous to the one in the Joseph- 
son junction array. The ±1 chirality spin can be immedi- 
ately read from the value of the (counterclockwise) phase 
change around each triangle ±7r, i.e. from the value of 
the induced flux through each triangle: or 1 flux quan- 
tum. Even though this realization seems quite similar to 
the previous one, there are differences that arise mainly 
from the fact that time-reversal is explicitly broken by 
the external field in the wire networks. For example, the 
±7r chiralities do not have the same energy in the case of 
wires of finite width. We refer the reader to the thorough 
discussion of the energetics by Park and Husei^. 

D. Mapping to a constrained Ising model 

The hard constraint of the 3-color model imposes a 
hard constraint in the allowed configurations of the chi- 
rality ±1 Ising spins. Here we show that in the spin 
representation the hard constraint requires that any ele- 
mentary hexagonal plaquette P must have total magne- 
tization: 

°p = J2 a * = ±6 > °- w 

iGP 

A similar result was obtained by Di Francesco and Guit- 
ter when connecting the folding problem in the triangular 
lattice to the 3-coloring modeU^. In our proof, we make 
use of phases accumulated along paths on the hexago- 
nal lattice, requiring that these phases are single valued. 
This approach is more appropriate to the discussion of 
superconducting systems and their defects (integer and 
fractional vortices) that we present in this paper. 

Indeed, as we show, one can obtain a simple interpre- 
tation of the hard constraint by identifying the accumu- 
lated phase around any loop lying on links of the hexago- 
nal lattice with the circulation of a vector potential. For 
concreteness, we will use the example of the Josephson 
junction array in the discussion, but the argument is gen- 
eral. 

The phase 9i ia on the edge a of the superconducting 
triangle i can be written as: 

Qi,a — @i + ei, a ' ^i,oi (2) 

where ei >0 is the unit vector that points from the center 
of triangle i to its a-th edge, and the "gauge" potential 
A i: a is defined at the center of such segment (see Fig. 121. 

The phase difference across a face a between triangles 
i and j is: 

di. a ~ dj.a = Qi — 0j + [h,a 1 M,a ~ &j,a 1 -4j,a] ■ (3) 
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The last term is simply the discrete sum equivalent of 
J d£ ■ A (notice that for neighboring sites i,j the unit 
vectors are opposed, e^ a = —e.j, a )- 

Now recall that one can write Qi^ a = 9i + -j-acj and 
hence the vector potential is such that: 

r 2n 

et, a • Ai, a = "y afJ i (4) 

What is the corresponding magnetic field? This is 
simple to answer, by looking at the accumulated phase 
around a loop. Consider an elementary anti-clockwise 
hexagonal loop. The loop visits six triangles, and the 
portion of the loop within each triangle enters through 
face a and exits through face a — 1 (mod 3), so that the 
accumulation of the vector potential along that portion 
of the loop is: 

ei,a-i • -Ai jQ _i - e i>a ■ A l:a = — (a - 1) <r» ^-acrj 



The above result, that each of the six sites visited by an 
elementary hexagon loop contributes — ^r-<Xj to an anti- 
clockwise accumulation of phase around the loop, has 
a very simple interpretation. Each Ising spin er s ; = ±1 
corresponds to a =p27r vortex sitting at a vertex of the 
hexagonal lattice. Each vertex is shared by 3 hexagons; 
hence each hexagon can be thought to contain 1 /3 of that 
vortex, as depicted in Fig.0J This is why the contribution 
from the hexagonal path going through vertex i picks up 
the phase — ^-<Ji as shown above. Basically, the vortex 
is divided equally among the three neighboring hexagons 
sharing the common vertex. 



sites. Now, matching the color scheme after going around 
any closed loop requires the phase around any hexagon to 
be uniquely defined (mod 2ir) , which in turn requires the 
flux to be a multiple of 2ir: 2-k/2> <jp = (mod 2ir), that is 
dp = ±6, (notice that dp is even). Since the total flux 
inside any loop is given by the sum of the fluxes through 
each elementary hexagon, then the condition <jp = ±6, 
grants the the phase to be uniquely defined (mod 2tt) 
around any loop. 

Once the dp = ±6, constraint is satisfied, there is 
a one-to-three mapping of any spin configuration to a 
configuration of the color model, since there are three 
even permutations of the colors that produce the same 
chirality spin configuration. 

In the case of the kagome wire networks at half-flux per 
triangle (or vertex of the hexagonal lattice), each trian- 
gle will accommodate either or 1 vortex. So instead of 
<Tj = ±1 one has a variable n% — 0, 1. Still, the vortices 
are split equally into three pieces, and the circulation 
around a hexagonal plaquette P going through the cen- 
ters of the kagome triangles is ^-Np = ^2 ieP n,i. The 
circulation is a multiple of 2ir if = 6, 3, 0. Indeed, 
the fact that the vortices in the elementary triangles are 
shared by three sites was used by Park and Husei^ in 
their argument for fractionalized vortices in the kagome 
superconducting wire networks. 

For finite U, there are defects that violate the ap = 
±6,0 constraint; we shall discuss these defects in detail 
in section llVl where we study integer and fractional vor- 
tices, as well as open segments of closed two-color loops. 
We analyze whether these different defects are confined 
or deconfined, and their importance in determining the 
ilk of the processes responsible for the dynamics. 




E. Interactions 

Each experimental realization of our model contains 
sub-dominant effects that may lead to a degeneracy lift- 
ing of the ground state. In this paper we concentrate 
on the effect produced by nearest-neighbor interactions 
between the chirality spins: 



FIG. 4: A vortex sitting at each vertex in the hexagonal 
lattice is shared by three hexagons. Hence, the contribution 
to an anti-clockwise accumulation of phase around a hexagon 
encloses one third of each of the six vortices sitting at the six 
vertices in the loop. 

Using equation JSJ we can now compute the flux en- 
circled by an elementary hexagon on plaquette P; it is 
given by: 

= -2tt/3 ^ <?i = -2?r/3 erg (6) 

Therefore the flux enclosed by an elementary hexagonal 
loop is just 1/3 of the sum of the vorticities in the six 



H=-^Ja i a j , (7) 

where the coupling J depends on the microscopic details 
of the problem. Such a coupling can arise, for example, 
if one considers the higher order effects of having an ex- 
tended Josephson junction barrier between two neighbor- 
ing triangles in the array geometry. In the Appendix [S] 
we show how to derive the constants U and J from a 
microscopic Hamiltonian for the array of Josephson cou- 
plings and we discuss the conditions for having U >• J. 
The sign of the J coupling is positive in this case. 

This nearest-neighbor interaction leads, in the color 
language, to an aligning or anti-aligning interaction be- 
tween the bonds, depending on the sign of the coupling 



6 



constant J as it can be easily seen with the help of Fig. [2] 
For J positive, the spin interaction is ferromagnetic and 
the zero-temperature ground state (GS) of the system has 
all the bonds with the same color aligned in the same di- 
rection. We will refer to this translation invariant state as 
the FMFS state, or single crystal state. For J negative, 
the spin interaction is antiferromagnetic and the zero- 
temperature GS of the system is a configuration where 
the six bonds in every hexagon form a sequence of only 
two alternating colors, which is simply the Neel order in 
the hexagonal lattice. 

In the following section, we discuss the thermodynam- 
ics of this system considering only the phase space of the 
configurations allowed by the ABC coloring constraint 
or, equivalently, by the Cp = ±6,0 constraint. 



III. THERMODYNAMICS OF THE 
DEFECT-FREE MODEL 

A. Small J and the CFT description 

Since the model without interactions can be described 
by a WZNW CFT, it is tempting to use this technique 
to analyze its behavior for small values of the spin-spin 
interaction. 

The first step is to represent the system by a height 
model (see Kondev et al. for details 5 ). Flat configu- 
rations of this height model correspond to the different 
Neel states of the system. In terms of the colors there is 
a total of six of those configurations which are arranged 
to form an hexagonal lattice. The coarse grained version 
is described by two fields h — (hi, /12) and a locking po- 
tential V(h) that favors the fields to lie in one of the flat 
configurations; this potential has then the periodicity of 
the hexagonal lattice. The action reads: 



S = 



d 2 x 



(Si 



V/i| 



V(h)) 



(8) 



In this language, the spin-spin interaction introduces a 
perturbation which is proportional to the "locking po- 
tential" since, depending on the sign of J, it favors or 
opposes the locking in one of the flat configurations. In 
the WZNW language, the locking potential can be writ- 
ten as a current-current perturbation of the underlying 
WZNW model 5 .. 

When the spin-spin interaction is turned on, we can 
use this description to propose an action for the per- 
turbed CFT. Since the A,B,C permutation symmetry is 
preserved, we can argue that the perturbing term to the 
pure CFT action should read: 



Mj2 J R J L a3 + J R ajJ L)\ (9) 



where the ctj's are the generators of the root lattice of 
su(3), and the Cartan generators J Hi are simply given by 
the derivatives of the height fields dhi. The case Xe = Xh 
corresponds to the su(3) symmetric case. The one loop 
Renormalization Group (RG) equation in this case reads: 



2tt 



(10) 



and for A > the flow is toward the unperturbed level 
1 su(3) WZNW model, which can be identified with the 
J = case. In general, however, we just have the A,B,C 
permutation symmetry, and we cannot exclude the pos- 
sibility of Ah 7^ Xe- Defining <5A = A# — Xe, the RG is 
now: 



SX 
X E 



- SX X E 

-i-A| -iw A 

Z7T 7T 



(11) 



where, at least for a small spin-spin interaction, we as- 
sume 1 5 X I << Xe- The RG flow is as follows (see Fig. [SJ. 
For SX > 0, the system flows to the line of fixed points 
Xe = 0. While the su(3) symmetry is broken, the system 
remains critical. We propose that this case corresponds 
to a ferromagnetic interaction, since it is equivalent to 
a decrease of the locking potential. This result is valid 
for small inter-spin couplings. As we show below, for 
large enough couplings a first order phase transition takes 
place. Since this is highly non-perturbative in the CFT 
language, this scenario is much better described by the 
Cluster Mean Field Method that we explain below. For 
an antiferromagnetic coupling, SX < and the flow goes 
toward strong coupling, bringing the system off criticality 
and forcing the system into anti-ferromagnetic ordering, 
as was argued by Huse and Rutenberg»s in their studies 
of the related classical kagome XY model. 



B. The Cluster Mean Field Method: general 
approach 

The CMFM is a technique that has proven to be very 
powerful in studying structural phase transitions in crys- 
tals and the thermodynamics of Vertex models^. When 
a system is constrained, fluctuations are considerably re- 
duced and an appropriate mean field treatment can give 
very good results if the constraint is taken into account. 
The idea is to consider as the fundamental entity cou- 
pled to a "molecular" field, instead of a single spin, a 
cluster in which the allowed spin configurations are re- 
stricted by the constraint. The bigger the cluster, the 
more accurately fluctuations and constraints are taken 
into account. This method has given very precise results 
for the ice model^ and is a good candidate for giving 
an accurate picture of our constrained spin model in the 
hexagonal lattice. 
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FIG. 5: Diagram of the RG flow for our model, where the 
horizontal axis corresponds to 5X and the vertical to Xe- The 
solid lines are numerical solutions of the system of equations 
llll l for three different initial conditions, and are drawn for 
visualization purposes only. 



It is particularly simple to introduce the CMFM in the 
case of a corner sharing plaquette*^ lattice with Hamil- 
tonian: 



• /i£<?t 



(12) 



ho 



where the range of the Jjj- interaction is shorter than the 
distance between the two farthest spins in a plaquette. 
This is the case for the present system. Let us assume 
that the lattice has N spins and 2N/S plaquettes, where 
each plaquette has S sites. The sums in the Hamiltonian 
can be rearranged as: 



i,jeP 



ieP 



hJ2^, (13) 



where the first sum is over all plaquettes P and the last 
term compensates for the double-counting of the site en- 
ergy term. The mean field approximation is obtained 
by considering each term as the sum over an elementary 
cluster (of S and 1 spins respectively) coupled to an ef- 
fective field representing the interaction with the rest of 
the lattice: 



H 



2N 
~S~ 



2N 



£ JijViVj + (h + 4>ext) £ Oi 

i.jeP ieP 

-N[(h + ^)a t ] 

-Hs-NHu (14) 



where H$ and Hi are the S- and 1-spin cluster Hamil- 
tonian respectively. Here, <j) and (j) ext are proportional to 
the number of spins that are external to the cluster but 
connected to the internal spins. Since for the 1-spin clus- 
ters such number of external spins is twice the number 
for the iS-spin clusters, we have <j) = ^4>ext- Let us now 
define the effective internal energy per spin: 



e=-{Hs)s-{Hx)x, 



(15) 



where (. . . }s an d (■ ■ ■ )i are the thermal averages com- 
puted with H$ and Hi respectively. Integrating then 
over the inverse temperature (3 we get an effective free 
energy: 



(3F: 



InZi, 



(16) 



where Zi = Tr{exp(— (3Hi)}, i = S,X, and the integration 
constant has been chosen such that in the case of uncon- 
strained spins we get the trivial entropy ln(2) at infinite 
temperature. Minimizing the effective free energy with 
respect to 4>: 



dF 



= 



(17) 



is equivalent to imposing the self-consistency equation for 
the magnetization: 



(<r)s = (o)] 



(18) 



and it gives us the optimal value for the field 0, which 
determines the behavior of the system at a given tem- 
perature. An important benefit of this method is the 
fact that it can be extended to larger and larger clusters. 
This allows to improve systematically the accuracy of the 
results. 



Application of the CMFM to the defects-free 
model 



In order to be able to apply the CMFM to our prob- 
lem in a straightforward way, it is convenient to switch 
to a bi-dual representation and describe our system in 
terms of spins Sij — ±1 sitting on the links of the hexag- 
onal lattice (see Figure 01. These spins are given by the 
product of the original chirality spins <Xj at the two ver- 
tices of each link: Sij — Oi<jj. Obviously, the number 
of configurations of the S spins is half the number of 
original a spin configurations, since we have quotiented 
by the global Z2 symmetry of the original model. The 
advantage of this mapping is that our lattice becomes 
now the (corner sharing hexagons) kagome net in which 
each spin Si is shared by two elementary plaquettes. In 
this description, the Hamiltonian (J7J restricted to the 
nearest-neighbor interaction reads simply: 



H 



(19) 
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where the index a refers to a link of the hexagonal lat- 
tice or a site of the bi-dual kagome lattice. The CMFM 
implementation is particularly easy since in this picture 
we just have an effective magnetic field J in (|f 2fl . The 
clusters that we use are the single spin cluster and the 
elementary hexagon cluster (with f 1 different configura- 
tions for the S spins), and the corresponding partition 
functions are given by: 

Z\ = ax 2 + l/(ax 2 ) 

Z 6 = a 6 x 6 +a~ 6 x- 6 + 3(ax) 2 + 6/(ax) 2 , (20) 

where a = e^ J and x = e^/ 2 . We can now obtain the 
values 4> pt corresponding to the minima of the effective 
free energy. Notice that (j) op t determines the equilibrium 
value of (S), i.e. of the internal energy per link of the 
original system. This method predicts the following sce- 
nario: for T — > oo we have (S) = f/3, which corre- 
sponds to an antiferromagnetic coupling in the system 
solely due to the constraint. This non-trivial value of the 
energy density is very close to the result obtained with 
the numerical method (see Sec. 0. The cluster mean 
field method also gives a reasonable estimate for Bax- 
ter's entropy in the limit T — » oo. Replacing 12Ufl in 
and taking the limit T — > oo we obtain the entropy per 
site S = ln(ll/8)/2 ~ 1.1726, while the exact value is 
1.2087.... Since the analytical expressions for the forth- 
coming quantities are too cumbersome, we just mention 
here their numerical values. At T ~ 9.872J the system 
undergoes a first order phase transition in which the en- 
ergy density jumps from (S) ~ 0.05 to a fully polarized 
state in which (S) is exactly —1 (see Figure|5J|. This tran- 
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FIG. 6: Plot of the internal energy per link as a function 
of the temperature. The solid line is the prediction from the 
CMFM and the doted-dashed line is the result from the nu- 
merical simulation. 

sition has been first noticed via transfer matrix analysis 
by Di Francesco and Guittei^ in the context of a folding 
transition. Our CMFM result is very close to their esti- 
mated critical temperature (9.1J) and even closer to our 
Monte-Carlo estimate 9. 6 J (see Sec. IV A"TJ) . In terms of 



the original spins, this behavior corresponds to the ex- 
otic scenario in which the magnetization jumps from to 
the fully saturated value 1 at the critical point, as was 
argued by Di Francesco and Guitter— . A similar kind 
of transition is also found in a frustrated spin model on 
the triangular lattice 1 ?, which turns out to be equivalent 
to a dimer model on the hexagonal lattice. Such kind of 
transition is accompanied by slow dynamics and aging. 
As we will see below, slow dynamics is also a central issue 
in our case. 

Another temperature that we can compute via the 
CMFM is the spinodal temperature of the system. This 
is typical of first order phase transitions, where an appro- 
priate fast cooling process can avoid crystallization and 
bring the system into a supercooled liquid phase. The 
spinodal temperature T sp is the temperature at which the 
supercooled liquid becomes unstable due to the crystal 
nucleation process. In this case, we can study the shape 
of the CMFM effective free energy as a function of </> for 
different temperatures. Starting from T ~ oo and lower- 
ing the temperature, the minimum corresponding to the 
liquid phase first becomes a local minimum (metastabil- 
ity) and eventually disappears. This metastability limit 
corresponds to the spinodal temperature T sp ~ 7.56 of 
the present model. 

The choice of the bi-dual spin representation to im- 
plement the CMFM is due to the fact that the system 
becomes a model for which the CMFM is particularly 
suitable. Indeed, in terms of the bi-dual spins, the sys- 
tem becomes a kagome lattice seen as an array of corner 
sharing hexagons, in which now the new spins are sitting 
at the vertices. By associating to each of the 11 config- 
urations for each hexagon its corresponding energy, the 
model can also be described as an 11 vertex model on 
the triangular lattice dual to the hexagonal. This choice 
of variable usually limits the analysis since it does not 
allow to measure the magnetization of the system, which 
is the typical order parameter used to study phase tran- 
sitions. In the present case however the energy density 
variable gives very good results in the characterization 
of the system since the transition is first order. For con- 
tinuous phase transitions the situation is different. Even 
though the CMFM still gives a quite accurate result for 
the numerical value of the energy density (in contrast 
to the normal mean field method), it may fail in repro- 
ducing a subtle behavior such as an infinite slope point 
at T c in the energy vs temperature curve. In this case, 
measuring the magnetization of the system is a much 
more powerful tool to detect and study the second or- 
der phase transition. Thus, one needs to get back to the 
original spins instead of the bi-dual ones. Implementing 
the CMFM technique within the context of the real spins 
has two main disadvantages in our case. On one hand, 
the spins do not form corner sharing plaquettes, and re- 
lating the mean fields acting on the 1-spin cluster and on 
the 6-spin cluster becomes more difficult. On the other 
hand, since the coupling J is a two spin nearest-neighbor 
interaction, a single variational mean field can not take 
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simultaneously into account both the J interaction (for 
which each spin interacts with its three neighbors) and 
the effective interaction due to the constraint (for which 
each spin interacts with all the 12 spins belonging to the 
three adjacent hexagons). 

D. Free energy argument for a first-order phase 
transition 

The key point for understanding this particular phase 
transition is to understand the very peculiar nature of 
its FMFS ground state. As we already discussed before, 
in the FMFS state all the bonds of the same color are 
aligned in the same direction. As a result, any two-color 
loop is maximally straight and winds around the whole 
system. Thus, the smallest possible rearrangement of 
the FMFS configuration that produces another allowed 
configuration is the update of one of such loops. This 
is a striking feature of the ferromagnetic three-coloring 
model: the GS is separated from the first (1-loop) "ex- 
cited" state by a system-spanning update which costs an 
energy: Ei-i oop - £fws = 2JL, where -Efmfs = -3JL 2 
and L is the system size (21? sites, 3L 2 bonds). Notice 
that if one prepares the system in the T = FMFS and 
starts to heat, the system is likely to remain in that state 
even for T — > oo for fast enough heating rates. Indeed, 
such an energy separation is likely to make the FMFS 
state metastable even for T — > oo, in the thermodynamic 
limit. Since the FMFS state has zero entropy and the en- 
tropy of a straight winding loop is ln(3i), we can write 
the free energies of the two states: 

-Ffmfs = — 3JL 2 

fkoop = -3JL 2 + 2JL-T -ln(3L). (21) 

Clearly in the thermodynamic limit the energy cost 
AE r~j L overwhelms the entropic gain AS ~ lni and 
the excited state will never be favored over the FMFS 
state at any temperature. A similar argument applies to 
higher excited states, as long as their entropy is not ex- 
ponential in the system size. The system is incapable (at 
equilibrium) to move out of its ground state in a "smooth 
way" . In terms of configurations, it has to jump from a 
fully ordered state into a state with finite domain size. 
Since it is reasonable to assume that a finite-domain-size 
configuration has negligible magnetization, we can intu- 
itively understand the origin of the complete first-order 
phase transition observed with the CMFM. 

The peculiarity of this transition and the relatively 
small variation of the internal energy in the disordered 
phase make it possible to obtain an estimate for the tran- 
sition temperature by comparing the free energy of the 
FMFS configuration with the free energy of the disor- 
dered configuration. In order to compute the free energy 
of a disordered configuration, we use the average infinite- 
temperature internal energy of the system Eoo = JL 2 , 
an estimate derived via the CMFM in the previous sec- 
tion and confirmed by the numerical results (see Sec. lVjl . 



Then, we can use Baxter's exact result for the residual 
entropy as an estimate of the entropy and obtain the free 
energy of a disordered state at all temperatures: 

F disordorod = J I? - T ■ 2L 2 In (1.2087). (22) 

By comparing the free energy of the FMFS state 
■Ffmfs = — 3JL 2 with -^disordered we obtain an estimate 
for the transition temperature 2 J/ (In 1.2087) ~ 10.55 J, 
which is reasonably close to the result from the CMFM 
T c ~ 9.872J. 



IV. DEFECTS AND THEIR ROLE IN THE 
DYNAMICS 

In this section we discuss the importance of defects in 
determining how the system can, dynamically, move from 
one of the allowed low energy configurations to another. 
For concreteness, let us start by discussing the Josephson 
junction arrays, i.e. the case of TL-i x U(l) symmetry. 

A. Integer vortices 

For finite U, it is best to understand the system in 
terms of the chirality Ising spins, plus XY spin waves of 
the U(l) sector. The lowest energy excitations over any 
configuration with Ising spins satisfying <7p = ±6, are 
topologically trivial (no vortices) XY spin waves. 

When the cxp = ±6, is preserved, vortices of the U(l) 
sector can only have vorticity that is an integer mul- 
tiple of 2n. These vortices cost an energy of order of 
magnitude U, the vortex core energy. The U(l) phase 
twist leads to the usual logarithmic interaction between 
a vortex/anti- vortex pair, 

£ 1 (xU2tt\uR. (23) 

and these pairs are confined below a Kosterlitz-Thouless 
type transition at a temperature scale T^j, oc U. Since 
we are interested in the regime of temperatures T <C 
U such that the three-color constraint is enforced, these 
integer vortices will be confined. 

Now, what are the accessible excitations that break the 
<7p = ±6, constraint? 

B. Fractional vortices 

A fractional vortex excitation is illustrated in Fig. [7| 
Such fractional vortices are always created in pairs via a 
nearest-neighbor exchange of opposite pointing spins and 
they have been discussed by Park and Husei^ in the case 
of the superconducting kagome network. A fractional 
vortex excitation corresponds to a single hexagon that 
violates the <7p = ±6, constraint. We define its frac- 
tional vorticity as T — 2n v — ap (mod 2ir). Thus, 
we have v = ±1/3 for Up = ^2 or <Tp = ±4. 
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FIG. 7: A pair of ±1/3- vortices created by a nearest-neighbor 
spin exchange. The solid-line lattice represents the kagome 
network considered by Park and Huse--. The Josephson junc- 
tion triangular array is represented instead by the bold tri- 
angles. The corresponding hexagonal lattice in our model is 
shown only around the two defects (dashed line). In the bot- 
tom right part of the picture we show the mapping to the 
bi-dual representation used in the CMFM. 

The presence of defects causes a fractional accumu- 
lation of the link sum of the vector potential Ai >a , the 
equivalent of <f> dl ■ A in the continuum limit, that equals 
±27r/3 (mod 27i"). Once again it is useful to resort to the 
picture in Fig. ^ to understand that only one third of 
the vorticity associated to an Ising spin at a vertex is in- 
cluded in the circulation around an elementary hexagon, 
and hence the flux is | • 27r o-p. 

To minimize the energy cost across the Josephson junc- 
tions, the superconducting phases 9i in the triangles 
must adjust accordingly to pick this extra phase differ- 
ence ±27r/3. Hence, an excited state that breaks the 
(Tp = ±6, constraint in the Ising sector must be ac- 
companied by a t/(l) phase twist that scales with the 
distance r from the defect as l/(3r) (in units of the lat- 
tice spacing). 

The U(l) phase twist leads to a logarithmic interaction 
between a fractional vortex/anti- vortex pair a distance R 
apart: 

£ 1/3 <xU^]nR. (24) 

Thermodynamically, there is an entropic contribution to 
the free energy, which was calculated by Moore and Lee^, 
and shown to also be logarithmic. Therefore, there is a 
confining transition of the Kosterlitz-Thouless type at a 
temperature U/9. If the Josephson coupling 

U is large compared to the temperature T, which is the 
regime we are interested in, then one is deep in the con- 
fined phase, and fractional vortices are rather ineffective 
as a source of phase space reconfigurations. 



C. Open segments of closed two-color loops 

There is a special way to flip Ising spins along certain 
strings lying on the hexagonal lattice that, while violat- 
ing the dp = ±6,0 constraint, only costs energy at the 
extremities of the string, irrespective of its length. 

To understand these excitations, let us start by looking 
at the simple case of a single spin flip that violates the 
constraint on three neighboring hexagons. In terms of 
the color model, all colors remain perfectly well defined, 
with the exception of the one vertex where the spin flip 
occurred. The energy cost of this defect is of order U. 
It is possible that locally adjusting the U(l) phase near 
the defect might slightly relieve this cost, but we have 
not investigated this issue. A single spin flip could split 
into a ±1/3 and ±2/3 (or equivalently, a —1/3) fractional 
vortex pair. These, however, are confined together at low 
temperatures compared to U, as we argued above. 

In the three-coloring model, this spin flip defect cor- 
responds to the initial step of creating an open segment 
defect described hereafter. Out of the three bonds de- 
parting from the spin flipped site, two must have ex- 
changed color (in order to change the chirality of the 
vertex), thus violating the color matching with the cor- 
responding two neighboring sites. If we now move these 
two color defects starting from the two neighboring sites 
and performing the same original color exchange, we can 
propagate the defects at zero energy cost along a pre- 
defined path. Indeed, every color exchange will fix the 
previous color mismatch and create a new one, one lat- 
tice spacing apart. Notice that this process will flip all 
the spins between the two end points along the path. It 
is useful to recall the color description of the allowed low 
energy states. Imagine one follows an ABAB... sequence, 
that always forms a closed loop in an allowed configura- 
tion. We have already seen that flipping the whole loop 
to BAB A... maintains the system in an allowed config- 
uration. It is also trivial to show that this update flips 
all Ising spins visited by the loop. While this is a rather 
non-local move, starting from a single spin flip (color ex- 
change) and propagating the color defects as above, we 
can realize this move through a sequence of local updates. 
Instead of flipping the whole loop at once, one can do it 
in steps, flipping the spins along a piece of the loop se- 
quentially. Notice that the energy cost of this string is 
paid only at the end-points and is of order U, as long as 
the sequence of spin flips moves on its two-color track. 
The end-points can be thought of as a defect pair con- 
nected by a string. This special path is hidden in the 
constrained Ising representation, but is clear in the 3- 
color one (see Fig. |SJ. The defect pair, once formed, 
can diffuse around the one-dimensional loop, and it has 
two channels to decay back into an allowed state: cither 
the defects recombine by going around the whole loop, 
leading to the BAB A... configuration, or they recombine 
without winding around the loop back to the original 
ABAB... configuration. These are the defects considered 
by Kondev et al~. In the CFT description, they corre- 
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FIG. 8: Defect pair at the end-points of an open string, 
with end-points highlighted (shaded circles) and the relative 
two-color path (bold links) shown in a configuration of the 3- 
color representation of the model. The end-points can travel 
freely along the path via nearest-neighbor color exchanges, 
such as the one outlined by the double arrow. Eventually, the 
two end-points recombine by either exchanging all the bonds 
along the path, or by leaving them all unchanged. 

spond to vertex operators with conformal dimension 1/2. 
While, as we mentioned, for a fixed configuration of colors 
there is no confining force between pairs, an effective in- 
teraction appears because of entropic reasons, producing 
an algebraic decay with the separation distance for the 
partition function in the presence of such defects. How- 
ever, for the dynamics one is really interested in the cost 
for a given configuration. Therefore, the formation and 
recombination of these defect pairs constitute the main 
mechanism responsible for the dynamical evolution of the 
system. 

The defect formation time just enters as an over- 
all rescaling of the time steps for loop updates. Also, 
since the time it takes for the defects to move diffusively 
around the 1-D loop is algebraic in the loop length (and 
not exponential) we can neglect this correction and sim- 
ply treat the whole loop update as a non-local elementary 
move, now with a justified local origin. 



V. DYNAMICS 

In order to study the dynamic properties of the sys- 
tem, we use Monte Carlo (MC) simulation techniques of 
an N = 2L 2 -site hexagonal lattice (3L 2 bonds) with pe- 
riodic boundary conditions. As we discussed in Sec. [H] 
the choice of the single-step update is non-trivial due to 
the color constraint. In Sec. II VI we argued that the open 
segments of closed two-color loops are the main actors 
in the dynamical evolution of the system, based on en- 
ergy and confinement considerations. Thus, without loss 
of generality, we consider only loop updates as single- 
step updates of our MC technique. We also assume that 
the rate of formation of the open segment defects is low 
enough not to allow for defect proliferation (i.e. for the 



intersection of two different open segments before they 
recombine) . 

To implement a loop update we proceed as follows: we 
first choose one site and two colors at random; then we 
compute the energy difference in the system for the up- 
date of the corresponding loop; eventually we accept or 
reject the update based on the usual Boltzmann proba- 
bility. Notice that, with this choice of the single MC step, 
the update of a loop takes one unit of time, independent 
of its length. In a possible experimental realization we 
expect the two ends of an open segment defect to walk 
randomly along the corresponding closed path, until they 
recombine. Thus, our MC dynamics is accelerated and 
the rescaling of our MC time with respect to a possible 
"real" time is highly non-trivial. Since we are interested 
in studying the slowing down and freezing of the dynam- 
ics in the three coloring model, we choose to use the 
accelerated loop dynamics in order to be able to sample 
much longer time scales, otherwise inaccessible with a re- 
alistic update mechanism based on defect formation and 
recombination. 

In terms of the loops, one can notice that the two or- 
dered configurations FMFS and Neel (ferromagnetic and 
antifcrromagnetic respectively) correspond to the two ex- 
trema in loop curvature. In the FMFS configuration, the 
loops are completely straight loops, winding around the 
whole system. In the Neel configuration, the loops are 
maximally curved into single hexagon loops. For these 
reasons, we expect an entropic jamming in the approach 
to the FMFS state, for a ferromagnetic choice (J > 0) of 
the interaction, as discussed in the case of infinite-range 
interactions by Chakraborty et al. A . Indeed, entropy fa- 
vors rough and entangled loops, which in the infinite tem- 
perature limit have a fractal dimension equal to 1.5^*i£. 
This creates a phase-space bottleneck due to the small 
number of configurations that allow the system to reach 
the FMFS state with straight, packed loops. On the other 
hand, the approach to the Neel state in the antifcrromag- 
netic interaction case ( J < 0) is much smoother for the 
system. Even though this state has zero entropy by itself, 
single-hexagon flips allow the system to achieve a gain in 
entropy of the order of In L 2 with an energy cost of the 
order of 6J. Indeed the Neel state corresponds to the 
ideal states defined by Kondev and Henley^, which have 
maximum entropy density in the sense that they allow for 
a maximum number of local rearrangements of the spins 
in accord with the constraint. Thus, we do not expect 
any jamming phenomena to play a role in this case. 

In this section we consider only the case of ferromag- 
netic interactions and we set J = 1 as the unit of mea- 
sure of energies and temperatures. In order to be able to 
access large simulation times, we choose the smallest sys- 
tem size for which our results do not show a significant 
dependence on system size {L = 18). 
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A. Transition temperatures 



We can then integrate to obtain the free energy: 



1. Estimate of the thermodynamic transition temperature 

The first result that we observe both in cooling/heating 
simulations and in quenching simulations is the phase 
space "isolation" of the single-crystal phase or FMFS. 
Even though at equilibrium the system must eventually 
favor the FMFS, we were unable to reach it within any 
simulation time, up to 10 7 MC steps. The system prefers 
to settle into a frozen polycrystalline (P-xtal) phase with 
zero or close to zero average magnetization, and with very 
slow, event-dominated dynamics. In Fig. we show the 
time evolution of the system after a quench in tempera- 
ture from T ~ oo to T = 6.0. After a single MC iteration 



(Fig. 9(a) |, only a few small crystalline seeds are visible 
in a disordered liquid background. These seeds quickly 
develop into well-defined domains (Fig. 9(b) I, whose size 
grows with time until the system b ecome s frozen into the 
polycrystalline (P-xtal) phase (Fig. 9(d) I. Notice the do- 
main boundaries following the "crystalline planes" of the 
hexagonal lattice in the polycrystal. The dependence of 
the crystalline mass m on time t reflects the remarkable 
slowing down in the dynamics once the system enters the 
polycrystalline phase. 

Even melting simulations starting from the FMFS 
phase and increasing the temperature are not useful to 
estimate the transition temperature. Indeed, they result 
in a large overestimate of T c , since the melting time re- 
mains much larger than the simulation time well above 
T c . 

The only measure we can achieve of the thermody- 
namic transition temperature T c is by computing the free 
energy in the liquid and crystal phases by integration of 
the internal energy. For a single crystal we know that 
/fmfs = — 1 at all temperatures, where / = F/(3L 2 ) is 
the free energy per bond. For the liquid phase, we use 
the curves in Fig. ll3l showing the dependence of the inter- 
nal energy on the temperature. Notice that the asymp- 
totic value of the internal energy at infinite temperature 
is different than zero. This is purely due to the con- 
straint, which appears to be slightly antifcrromagnetic in 
nature. A simple way to visualize this effect is to look 
at an infinite temperature configuration after performing 
a spin-flip operation on one of the two sublattices of the 
hexagonal lattice. The result is shown in Fig. I1UI 

An appropriate fit of the common high-temperature 
region of the internal energy (per bond) curves^: 



^liquid (T) 



a/T b 



(25) 



gives a ~ 4.3, b ~ 1.22 and c ~ 0.336. Notice that a 
naive high temperature expansion in powers of y may 
be plagued by the criticality at hight temperatures. In 
this sense the nontrivial exponent b may have an inter- 
pretation in terms of the CFT description at T — > oo. 



W) = /Wo) +[ d[3'£ (/?'); (26) 
J 13a 

setting /3q — for the liquid phase and using the known 
residual entropy of the system, we obtain: 



/liquid (T) 



ln(1.2087)T + c 



(b+ \)T b ' 



(27) 



where the 2/3 factor in front of the residual entropy 
comes from the fact that there are 3 bonds every 2 spins. 
Setting /liquid (T) = /fmfs = -1 gives the melting tem- 
perature T c = 9.6, in good agreement with the results 
from the CMF method. 

Even though T c is the actual thermodynamic transi- 
tion temperature, we are unable to observe this transi- 
tion due to the incredibly large time scales involved in 
the approach to the FMFS state. As it appears from the 
results below, the system seems to be completely unable 
to sample the phase space region corresponding to the 
crystalline phase, at least on our simulation time scales, 
and it is confined to an "effective phase space" . 



2. The dynamic freezing transition 

Instead of going through the thermodynamic transi- 
tion, the system remains in a supercooled liquid state be- 
low T c , until it reaches a temperature T* where it evolves 
into a frozen polycrystalline state. 

Looking at Fig. |9(f)| we can clearly see that the poly- 
cristallization is complete, in the sense that the domain 
boundaries are fully one-dimensional, with almost no in- 
terstitial liquid left. While the size of these domains in- 
creases with longer waiting times, the growth becomes 
extremely slow, basically stopped within our Monte Carlo 
time scales before reaching the single crystal configura- 
tion. This can be observed, for example, in the behavior 
of the zero-temperature saturation value of the energy 
in Fig. 13(b) and in Fig. The energy is in fact a 
measure of the area-to-perimeter ratio in the polycrys- 
talline phase, provided complete polycrystallization has 
been achieved. This is clearly the case in the T — > 
plateaus in Fig. Instead of approaching the value — 1, 
characteristic of the FMFS state, these plateaux seem to 
approach a limiting value £ F - xtaX (T = 0) ~ —0.74 for 
larger cooling times. 

The transition at T* can be seen as a dynamic phase 
transition and does not have a thermodynamic origin. 
However, we can reasonably establish a correspondence 
of this transition to a "true" thermodynamic phase tran- 
sition in a related, more constrained system. As we show 
with the following analysis, the origin of the dynamic 
transition at T* resides in a free energy barrier that 
prevents the system from visiting a phase space region 
around the FMFS phase, at least within our simulation 
timescales. Since only winding loop updates can change 
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(e) m = 0.68, t = 5.7 ■ 10 4 MC steps 



(f) m = 0.73, t = 5.4 • 10 5 MC steps 



FIG. 9: Time evolution snapshots of the system after a quench from T ~ oo to T = 6.0 (at time i = 0) below the transition 
temperature T* ~ 8.1. The dots represent the 2L 2 vertices of the hexagonal lattice (L = 36) and the two colors correspond 
to the two values of the chirality spin. The lattice is wrapped along the horizontal axis and along the 60° axis rotated 
counterclockwise above the horizontal. For each configuration, we report the measured crystalline mass m and the time t from 
the temperature quench. 
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(a) Original spins 



(b) Same spins after a one-sublattice spin-flip 
operation 



FIG. 10: A picture of a L — 18 system at infinite temperature: (a) the original chirality spins are shown; (b) we performed a 
spin-flip operation on one of the two sublattices in order to make visible the antiferromagnetic correlations due to the constraint. 



the number of bonds per color per direction, it is possi- 
ble to divide the phase space into topologically separated 
sectors by forbidding the update of winding loops. The 
FMFS configuration would then be in a topological sec- 
tor by itself, and starting from an infinite temperature 
configuration with equal number of bonds per color per 
direction it would be impossible for the system to reach 
its natural ground state. With this constraint, the sys- 
tem is expected to show a phase transition into a state 
which is not the FMFS, with a behavior analogous to the 
one observed in the present model. 

This polycrystal transition is an intrinsic transition of 
the supercooled liquid phase, which would not exist in 
the infinite time limit. If we were able to wait infinite 
simulation times, we expect the dynamic transition at 
T* to disappear, replaced by the equilibrium transition 
at T c > T*. 

Since we cannot apply the same technique used above 
for T c to the polycrystalline state, we have to measure 
T* with a somehow more empirical method. We first 
prepare the system into an almost completely polycrys- 
tallizcd state by cooling it at very low rates. We then 
chose a particular value for the temperature T and let it 
evolve in time. If it eventually reaches the liquid state, 
then we conclude that T > T*; conversely if it completes 
the polycrystallization process. The choice of the initial 
state closer to the polycrystalline state rather than to the 
liquid one is merely due to the stronger metastability of 
the liquid phase, as it appears from the asymmetry in the 



hysteretic process with respect to T* (see Fig. |l3(a) I . In 
Fig. ^] (top) we present the results in terms oitime evo- 
lution of the energy. Even though we do not have a sharp 
distinction between the behavior above and below T* we 



can clearly identify a transition at T* ~ 8.1 ± 0.1. When 
the system is set to a temperature T > 8.2, it quickly de- 
parts from the quasi-polycrystallyzed initial state, while 
for T < 8.0 it completes the polycrystallization process, 
thus lowering its energy. It is interesting to notice that 
all the quenching temperatures are below the thermody- 
namic transition temperature T c — 9.6, while the system 
behaves as if it is incapable of visiting the favored FMFS 
configuration. 



Since the total magnetization of the system remains 
close to zero for all temperatures and time scales that we 
are able to sample, it cannot be used as an order parame- 
ter for this transition. A more appropriate order param- 
eter is probably the crystalline mass to, shown in Fig. 1111 
(bottom). As proposed by Cavagna et alm&, the crys- 
talline mass measures the fraction of crystallized spins 
independently of the size of the polycrystals. We first 
define the elementary crystal unit as the four spin cluster 
composed by one spin and its three nearest-neighbors. To 
avoid double-counting, we choose the central spin exclu- 
sively in one of the two sublattices of the hexagonal lat- 
tice. Then, we define the crystal mass density m € [0, 1] 
as the number of these elementary units present in a given 
configuration, normalized by the total number of units 
L 2 . Since we need to keep the elementary unit small 
enough to be sensitive to small amounts of crystal mass, 
we have a limited power of resolution. In fact, even a 
random configuration has a non-zero average crystalline 
mass too — 0.01, which we consider as the effective zero 
of to. The results obtained by measuring the time evolu- 
tion of to are in good agreement with the conclusion that 
T* ~ 8.1 ± 0.1. 
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FIG. 11: Time evolution of the internal energy and crystalline 
mass, after the system has been prepared in an almost poly- 
crystallized configuration. The curves correspond to different 
quenching temperatures both above and below the transition 
temperature T* ~ 8.1 ± 0.1. Note that all the tempera- 
tures are below the thermodynamic transition temperature 
T c — 9.6, while the system behaves as if it is incapable of 
visiting the favored FMFS configuration. 



3. Some considerations on the dynamics of the polycrystal 

The data shown in Fig. are averages over 32 dif- 
ferent histories starting from the same initial configura- 
tion. The reasons for the large time fluctuations and the 
lack of a sharp distinction between above-T* and below- 
T* behavior, as shown instead in the system studied by 
Cavagna et alm&, are to be found in the peculiar, rare- 
event-dominated dynamics of the polycrystalline phase. 
It is worth to analyze this dynamics in detail, as it helps 
understanding also the phase-space isolation of the ther- 
modynamic GS, i.e. the FMFS crystal. 

With some simple reasoning about the colors and the 
chirality spins, one can see that within a single, ferro- 
magnetically ordered domain, all the bonds of the same 
color are aligned in the same direction. Thus, any two- 
color sequence inside the domain follows a straight path 



from one side to the other along one of the three crys- 
talline directions (or crystalline planes) of the hexagonal 
lattice. This high level of order is responsible for the 
first important difference with respect to usual domain 
growth: there are no small loops across the boundary of a 
domain (but for possible corner loops) and the domain is 
not capable of small rearrangements of its walls. While 
for example in a normal Ising model a domain can ex- 
pand gradually, in our constrained Ising model a domain 
can only crack from side to side. It is important to notice 
that these cracks will almost always bring the system into 
an excited state with higher energy, the energy difference 
being proportional to the length of the crack. 

If we now extend these considerations to the almost 
complete polycrystalline phase that the system is able to 
achieve below T* (see Fig. |?5J , we can see that any loop 
has to cross a few domains before closing on itself. In fact, 
bending of the loops are allowed only at domain bound- 
aries. Therefore, we have a second important difference 
with respect to usual domain growth: one domain cannot 
expand at the expenses of a single other domain; rather, 
the above cracks involve at least six domains (but for 
the case of winding loops), since every domain boundary 
corresponds to a 60° bending in the loop. One can easily 
convince oneself that the closer the system is to the poly- 
crystalline phase, the more the dynamics become frozen, 
requiring entangled, multiple-domain cracking in order to 
move from one configuration to another. This behavior 
can be seen for example by looking at the behavior of the 
spin-spin autocorrelation function (see ean. (|29[l ). shown 
in Fig. E| For small values of t W) the s ystem is still in 
a rapidly-changing liquid phase (see Fig. 9(b) I, and the 
correlation function roughly follows the stretched expo- 
nential behavior with a very short relaxation time dis- 
cussed in Sec. IV B 21 As the system gets deeper into the 
polycrystalline phase for t w — 2 ■ 10 3 or even more for 
t w = 2 • 10 4 (see Fig. 9(d)| ), the behavior of the correla- 
tion function shows how the system now evolves mostly 
via rare events that are responsible of extended changes 
in the system configuration. Notice the Z2 symmetry 
of the system. When the dynamics become highly en- 
tangled in the polycrystalline phase (see Fig. 9(e) I, the 
number of allowed configurations drops dramatically and 
rearrangements that bring the system from one configu- 
ration to its mirror image play a significant role in the 
evolution of the system (Fig. 



12(c)! 



It is important to underline the large energy cost of 
these updates, which scale with the linear size £ of the 
domains. Indeed, we can interpret this energy difference 
as the activation energy Ea{£) for domain growth. Pro- 
cesses where the activation energy depends on £, or more 
generally where freezing involves a collective behavior de- 
pendent on £ belong to classes 3 and 4 for growth kinet- 
ics^. In the next paragraph, we will address this classi- 
fication in greater detail. 

Even if the system is able to overcome the activation 
energy barrier, the three-coloring constraint plays a new 
key-role in preventing the system from reaching a new 
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FIG. 12: Spin-spin autocorrelation function C(t w ,t) for a single MC simulation and four different values of t w = 20, 2 ■ 10 2 , 
2 • 10 3 and 2 • 10 4 MC steps. The temperature is quenched at t = from T = oo to T = 6, the same used in Fig- El At t m — 20, 
the system is still in a rapidly-changing liquid phase (see Fig. |9(b)| . As the system gets deeper into the poly crystalline phase at 
t m — 2 • 10 3 or even more at t w = 2 ■ 10 4 (see Fig. |9(d)} , the behavior of the correlation function becomes discontinuous, reflecting 
a rare event do minated dynamics where the system undergoes highly non-local rearrangements. Notice t he Z2 symmetry of 
the system (Fig. |12(cj) . When the dynamics become highly entangled in the poly crystalline phase (see Fig. |9(e)| , the number 
of allowed configurations drops dramatically and rearrangements that bring the system from one configuration to its mirror 
image play a significant role in the evolution of the system. 



configuration. Let us consider an excited state after one 
loop has been updated in the polycrystalline phase. The 
system has then three types of updates available: the 
trivial repair of the crack, with consequent lowering of the 
energy; an independent update, which requires to over- 
come a similar activation energy; and the peculiar loop 
updates that are adjacent to the open crack. Clearly, 
since a loop update corresponds to flipping all the spins 
along the loop, the latter update has a vanishing energy 
cost because the original crack crosses crystalline ordered 
domains. Thus, the system is able, via these adjacent 
loops, to expand or contract a crack with essentially equal 
probability. Indeed we expect this process to be similar 



in nature to a random walk, with two possible outcomes: 
the crack eventually contracts and closes on itself, or all 
the domains involved in the original crack get essentially 
flipped, with minimal structural change in the original 
configuration. Notice that the last update in this pro- 
cess is of the repair type, with the system getting back 
to a lower energy state. The time to complete this pro- 
cess is the lifetime Td of a crack in the system, while the 
formation time of a new crack is determined by the ac- 
tivation energy barrier Tf ~ exp[— @Ea(£)]- At low tem- 
peratures, Td is much shorter than Tf] the system freezes 
into a specific polycrystalline configuration and the dy- 
namics involve only rare events where entire domains are 
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flipped simultaneously. At temperatures close to T* in- 
stead, Td becomes comparable to tj and multiple cracks 
allow the system to deeply rearrange the domains. Notice 
however that it is still a rare-event dependent dynamics. 
In a typical process of configuration change, the system 
visits highly excited states with complete "melting" of 
extended areas of the polycrystal, before freezing again 
into a new polycrystalline configuration . These highly 
excited intermediary states easily become long lived due 
to the metastability of the liquid phase, which has instead 
very fast dynamics (see Fig. [Poland the results hereafter). 



B. One-time quantities 

1. Energy vs temperature and growth dynamics 

In order to get a better insight in the dynamics of the 
model, we study the behavior of the system through tem- 
perature hysteresis with different cooling/heating rates. 
We vary the temperature from T = 40, where the liquid 
phase is stable and equilibrates very easily, down to T = 
and up again to T = 40, with a constant rate given by 
r = AT/ At, At being the total time to go from T = 40 
to T = 0. During these simulations we measure all the 
relevant quantities in our system: the internal energy, 
the magnetization, the staggered magnetization, and the 
crystalline mass. Both magnetizations remain close to 
zero for any temperature and cooling/heating rate. The 
behavior of the internal energy is shown in Fig. ED f° r 
some of the cooling/heating rates that we consider. The 
behavior of the crystalline mass is in agreement with the 
internal energy and does not provide any additional in- 
formation. 

The hysteresis observed in the energy curves is typical 
of first order phase transitions. From Fig. 



13(a) 



see that the hysteresis gets narrower for smaller values 
of r, indicating a transition temperature that is consis- 
tent with our previous estimate T* ~ 8.1 ± 0.1 (that 
estimate is also confirmed by looking at the position of 
the peaks in the specific heat, measured from the energy 
fluctuations, for different cooling/heating rates). Notice 
the asymmetry of the hysteresis toward the liquid phase, 
particularly evident for large cooling/heating rates, due 
to the metastability of the liquid with respect to the poly- 
crystalline ph ase. Fo r large cooling rates, see for example 
r = 0.4 in Fig. |13(b)| the energy curves never cross below 
the extrapolatecTTiiq U id(T) curve (dashed line in the fig- 
ure). Thus^S, the system does not polycrystallize and it 
remains in a supercooled liquid phase with respect to the 
polycrystalline phase until T = (recall that the liquid is 
already supercooled with respect to the FMFS phase for 
T < 9.6). This is confirmed also by the absence of a peak 
in the specific heat curves. As the temperature is low- 
ered to zero, the curves reach a final value of the energy 
that decreases monotonically with smaller cooling rates. 
But for very large values of r (larger than 0.4), this final 
value of the energy is reached already at a finite temper- 
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(a) Internal energy vs temperature hysteresis for 
three different values of the cooling/heating rate: 

r = 0.04, 0.004 and 0.00004. The hysteretic 
behavior is typical of a first order phase transition 
and is in good agreement with our measure of T*. 
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(b) Internal energy vs temperature plot for 
different cooling rates r = 0.4, 0.04, 0.004, 0.0004 
and 0.00004. The dashed line is the extrapolated 
internal energy of the liquid phase (eqn. I25l ~). 
Notice that, for r = 0.4, the system stays in the 
liquid phase till T = 0, since the energy curve 
remains above the dashed line at any 
temperature™. 



FIG. 13: Internal energy vs temperature behavior for our sys- 
tem, in the temperature range T G (0,40). These curves are 
obtained from simulations where the temperature is changed 
at a constant cooling/heating rate. For large temperatures 
(T > 15), all the curves overlap and the system is at equilib- 
rium in the liquid phase. Notice that there is no sign of the 
thermodynamic transition at T c = 9.6, as the system goes 
smoothly into the supercooled liquid phase. 



ature and the curves show a plateau typical of frozen or 
very slow dynamics. While we expect this behavior when 
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the system enters the polycrystalline phase, we can no- 
tice that this plateau is also present for curves where the 
system remains in the super cooled liquid phase (e.g. see 
the curve for r = 0.4 in Fig. 13(b) I. A detailed analysis 



of this behavior is beyond the scope of the present paper 
and will be addressed in the future. 

The dependence of the T — value of the energy on 
the cooling rate reflects the type of domain growth in 
the system. In particular, when the system enters the 
polycrystalline phase where domain boundaries are one- 
dimensional, the energy difference £(T = 0) — £fmfs = 
£ (T = 0) + 1 is proportional to the inverse of the linear 
size of the domains^. In Fig. [21 we show the behavior 
of £(T = 0) — £fmfs as a function of r. 
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confirmed, it implies that the behavior of our system for 
r G [8 • 10~ 5 ,0.2] belongs to class 4 growth kinetics^. 
Both class 3 (corresponding to the case of m — 1) and 
class 4 kinetics are typical of processes that involve a £- 
dependent collective behavior in the frozen phase. As 
discussed above, we indeed expected the system to show 
this logarithmic behavior. 

Eventually, for r < 8 • 10~ 5 the energy saturates to a 
limiting value (J 1 — 0) ~ —0.74, in agreement with 

the entropic argument we provided before. The system 
behaves as if a whole region of phase space around the 
FMFS configuration is dynamically inaccessible due to a 
very large free energy barrier. 

To further confirm this peculiar free energy landscape, 
we use again the CMFM described in Sec. [nin Fr om the 
numerical results, we assume as a first-order approxima- 
tion that the dynamically excluded configurations corre- 
spond to system energies smaller than the limiting value 
£P-xtal( T = 0) 0.74. We then impose appropri- 
ate constraints on the variational parameter such that 
the only allowed energies in the CMFM are larger than 
£ p " xtal (T = 0). Under these constraints, the method 
predicts a first order phase transition at T* ~ 8.36, in 
good agreement with the numerical value T* ~ 8.1 ±0.1, 
considering the approximations underlying this CMFM 
result. 



2. Domain nucleation vs liquid relaxation 



FIG. 14: Semi-logarithmic plot of the plateau value of the 
internal energy with respect to the GS energy of the perfect 
crystal (<?fmfs = — 1), versus the cooling rate r = AT /At. 
Three distinct behaviors can be identified: a power law be- 
havior £ ~ r 011 f or r > 0.2, when the system remains in 



the liquid phase; a logarithmic behavior £ 



ln(l/V u 



for 



8-10 < r < 0.2; and a saturation plateau at £ 



0) 



-0.74 for r < 8 • 10" 



Here we study the equilibration time of the liquid phase 
in comparison to the nucleation time for the polycrys- 
talline phase. 

We measure the connected piece of the two-times au- 
tocorrelation function: 



C(t w , t) 



1 

21? 



(29) 



As long as the system remains in the liquid phase, 
i.e. the energy curves never cross below the extrapolated 
^liquid CO curve, the energy follows a power law depen- 
dence on r: £ — £fmfs ~ r 011 . This is typical of class 1 
growth kinetics, where freezing originates from local de- 
fects with activation energies independent of the domain 
size 

As we lower the cooling rate, we reach a threshold 
where the energy curves start crossing the extrapolated 
^liquid (T) curve and the system polycrystallizes. This 
threshold happens at r t h — 0.2 and £th — —0.39. Be- 
low this threshold, the behavior of the energy changes 
abruptly into a logarithmic form: 



£{T = 0) — f FMFS = 



1 



1 + A 



^ ( r ■ Ti ' 



(28) 



From a fit of the results we obtain m ~ 0.85, even though 
our numerical data do not have enough accuracy to ex- 
clude the case m — 1. If our measurement of m ^ 1 is 



where (...) indicates the average over initial configura- 
tions of different MC simulations. Notice that J2i &i(f) — 
for all values of t within our simulation time scale, thus 
the disconnected piece of the autocorrelation function 
vanishes. Since we are interested in the relaxation time 
of the liquid phase at equilibrium, we quench the system 
from infinite temperature down to the target tempera- 
ture T and we wait for it to equilibrate. The correlation 
function becomes time-translation invariant and depends 
only on the time difference t — t w . At equilibrium, we ad- 
equately fit C(t — t w ) with a stretched exponential, which 
is the expected equilibrium behavior in supercooled liq- 

uids2a : 



C(i)=cxp [-{t/rf] 



(30) 



From the fit we obtain the relaxation time r as a func- 
tion of the quenching temperature, as shown in Fig. 1151 
We can extend the measurement of r below T* because 
of the metastability of the liquid phase. The system is 
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FIG. 15: Liquid phase relaxation time r as a function of 
temperature, as measured from the stretched exponential fit 
of the autocorrelation function at equilibrium. The dashed 
line corresponds to a power law fit while the dotted line cor- 
responds to a Vogel-Fulcher-Tamman fit. Notice that there is 
no dynamic signature of the transition at temperature T* . 



able to equilibrate as a supercooled liquid well before the 
polycrystal transition takes place, at least for tempera- 
tures close enough to T* . Notice that there is no dynamic 
signature of the polycrystal transition at T* in the liq- 
uid relaxation time. The Kohhausch exponent j3 of the 
stretched exponential fit decreases with temperature, as 
for realistic models of liquids. In Fig. ^] we show the fit 
of the r data both with a power law: 



.4 



-, T c lq = 7.0 7 = 1.0, 



(T - Tc q )t 

and with a Vogel-Fulcher-Tamman (VFT) form: 



t = tq exp 



A 



T-T Q 



T n =4.4 A = 11.1. 



(31) 



(32) 



The results of these fits have to be considered with ex- 
treme care. Because of the accelerated non-local dynam- 
ics and because of the onset of polycrystal nucleation, the 
temperature range where we are able to measure the re- 
laxation time of the liquid phase allows for r to vary only 
over a narrow interval, from 0.05 to 0.5 MC steps. As a 
consequence, the values obtained for the fitting parame- 
ters lack in accuracy, since the fit spans a single decade of 
data. Moreover, a VFT behavior typically involves the 
large r limit of the r(T) curve, which is not accessible 
in the present system due to the rapid nucleation of the 
polycrystal. Indeed, our numerical data are the tail of a 
possible VFT behavior, and they suggest that a VFT be- 
havior may be observed in the liquid phase of this system 
if the polycrystallization process were to be avoided. 

Since the correlation function decays to zero in ap- 
proximately 20 t, we can take this value as the equilibra- 
tion time for the liquid phase at a given temperature^: 
T cq(T) — 20 t(T). 



Measuring the nucleation time of the polycrystalline 
phase in this system is instead more complicated. Due to 
the frozen nature of the polycrystalline phase, we cannot 
compute its free energy as a function of temperature as 
we did for the liquid phase (see Sec. IV ATI) . Thus, meth- 
ods such as the one in Cavagna et alm& are not applicable. 
More naively, we have to estimate r nuc i directly observing 
the time evolution of the system. In Fig. ^| we plot the 
energy dependence on time for quenches of the system 
from infinite temperature to the target temperature T. 
As we discussed above, the system polycrystallizes when 
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FIG. 16: Time evolution of the energy of the system, after 
a quench from infinite temperature down to a target temper- 
ature T = 5.5, 6.0, 6.5, 7.0, 7.5, 8.0 and 8.5. The horizontal 
line corresponds to the energy threshold for polycrystalliza- 
tion £ t h — —0.39, as identified above. 

the energy falls below a threshold value E t h — —0.39. 
Here we use this value in order to identify the onset of 
the polycrystallization process in the energy curves in 
Fig. El The time when the system starts developing a 
polycrystalline phase is indeed the nucleation time T nuc i 
we are interested in. We can see that r nuc i ~ 800 at 
T = 7.5 while it drops to r nuc i ~ 170 at T = 7.0. 

Comparing these results with the ones of Fig. 1151 pro- 



vided we perform the rescaling r e , 



20 t, we can see 



that the crossover r cq = r nuc i will happen at a temper- 
ature T sp close to Tp, where the liquid relaxation time 
shows a rapid growth. We can reasonably locate this 
crossover in the temperature range 7.0 < T sp < 7.5. This 
temperature is the spinodal temperature corresponding 
to the metastability limit of the liquid, when the liquid 
equilibration timescales become of the same order of the 
nucleation timescales and the liquid phase becomes un- 
stable. The system reaches this limit in a time t sp of the 
order of a few hundred MC steps. 



VI. CONCLUSIONS 

In this paper we have studied the very interesting prop- 
erties of a model for describing the behavior, both static 
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and dynamic, of different arrays of superconducting de- 
vices. Among the examples discussed, the main candi- 
date to see such a rich phenomenology is a Josephson 
junction array of triangular grains of superconductors 
with p x ±ip y order parameter. In the limit of very strong 
Josephson couplings, the system is equivalent to Baxter's 
three color model in the hexagonal lattice. This model 
can in turn be represented by an Ising model with a con- 
straint on the total magnetization for each hexagonal pla- 
quette, Up — ±6, 0. In this paper we have presented a 
proof of this mapping based on the condition of the singl- 
evaluedness of a superconducting order parameter. The 
Ising degrees of freedom correspond, in the Josephson ar- 
rays with p-wave islands, to the chirality of the p x ± ip y 
order parameter. 

Within the constrained dp = ±6, space, the sys- 
tem is critical at infinite temperature but orders at any 
finite temperature if antifcrromagnetic interactions be- 
tween the Ising spins are present. For ferromagnetic in- 
teractions, it remains critical until a very particular first 
order phase transition takes place, where the system or- 
ders completely. This behavior is due to the peculiar 
nature of the ordered state, which is isolated in phase 
space from any of its excitations by an energy of order 
the system size. 

For a finite Josephson coupling strength, defects are 
present in the system, and there are violations of the 
color and, consequently, dp = ±6, plaquette constraint. 
A particularly interesting kind of defect is a fractional 
vortex pair. Within the context of the Josephson ar- 
ray of p x ± ip y superconducting islands, not only there 
is a large energetic cost to create these excitations, but 
they are also confined at low temperatures by logarith- 
mic interactions. The other kind of interesting excitation 
is formed by flipping the spins along open segments of 
closed two-color loops. While there is also an energetic 
cost to create them, these defects can circulate on the 
lattice without further energetic cost, in contrast with 
the fractional vortices. Moreover, a new defect-free color 
configuration is obtained through the process of creation 
of a string of spin flip excitations, the propagation of the 
defect along the two-color loop, and the recombination of 
the ends of the string after closing the loop. This mech- 
anism is precisely the microscopic origin of the Monte 
Carlo dynamics that we implement in this paper. 

Because of the constraint, the dynamics of the system 
is very peculiar. While the existence of a super-cooled 
liquid phase is typical of first order transitions, for our 
constrained system we find a whole temperature range 
in which such super-cooled liquid is stable for extremely 
long time scales. Indeed, at all the time scales studied in 
this paper, the fully order phase can not be reached and 
the system orders into a polycrystalline phase in which 
the global Z2 symmetry is unbroken. The transition from 
the liquid state to the polycrystal takes place at a critical 
temperature T* , smaller than the static (avoided) critical 
temperature. This dynamical temperature T* has been 
obtained both by studying the time evolution of the sys- 



tem after preparing it in a polycrystalline configuration, 
and by quenching the system from the liquid phase. The 
values obtained for T* are in agreement with the naive 
estimate that we are able to obtain from the CMFM tech- 
nique. The numerical analysis of the nucleation time and 
the liquid phase relaxation time allows us to give an es- 
timate of the spinodal temperature of the liquid. 

The rich phenomenology of the dynamics of this sys- 
tem is also reflected in the dependence of the difference 
between the final internal energy reached by the system 
and that of the fully ordered state on the cooling rate. 
While for very fast cooling rates this dependence shows 
a typical power law scaling, the nucleation of the poly- 
crystalline phase produces a logarithmic behavior until 
a total arrest in the domain growth is reached, meaning 
probably another logarithmic growth but with a much 
longer time scale. The origin of this scenario is the fact 
that the energy barriers through which the system has 
to pass to reach states with larger clusters grow with the 
size of the clusters. This places our system as one of the 
rare cases without randomness in which the dynamics is 
of class 3 or 4 in the classification of Lai et al~. 

An important open problem concerns the possible 
mechanism to get out of the polycrystalline state. Prolif- 
eration of other (confined) type defects, such as fractional 
vortices, is a possible mechanism to help overcome the to- 
tally arrested dynamics in the polycrystalline phase. In 
this case, the large time scale dynamics could be governed 
by the energetic cost of making a rather rare event domi- 
nated proliferation and circulation of such (confined) de- 
fects. It is noteworthy that, in the polycrystalline phase, 
not only the fractional vortices are confined (logarith- 
mically, with a prefactor of order U), but also the ex- 
citations that we argued are responsible for the micro- 
scopic dynamics, the open segments of closed two-color 
loops. The confinement of the two-color segments is pro- 
portional to the string length (linear) inside any ferro- 
magnetically aligned domain, with a prefactor of order 
J. The example that we studied in this paper suggests 
an interesting scenario where defect confinement at the 
microscopic level is responsible for the slow dynamics and 
out-of-cquilibrium behavior of a macroscopic system. 
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APPENDIX A: MICROSCOPIC ORIGIN OF THE 
U AND J TERMS 

In this appendix we estimate the relative values of U 
and J in terms of some microscopies for the tunneling 
through a Josephson barrier. Consider two neighboring 
triangles as in Figure sharing a common edge labelled 
by a. The microscopic tunneling Hamiltonian from a 
triangle labelled 1 to a neighboring triangle labelled 2 
can be written as: 



and 



H 



k,q 



h.C. 



(Al) 



Using second order perturbation theory, we can estimate 
from this expression the Josephson coupling between the 
two superconductors in a standard way. The result is: 



k . q 

where we used t? -* = t* 



E A 



(A2) 



k,—q 



_ and E\ is the superconduct- 



ing energy gap. 

It is useful to define the angles 



and 



as those 



formed by the vectors k, q and the reference unit vector 



ei,o- 



side labelled by a (see Fig. [3] for the definition of these 
unit vectors). 

The order parameters can be written as: 



<4 c -?> = A l 



= A (T^+^H) (A3) 

= Ae'( w< '*), (A4) 

where A is the order parameter magnitude, 6\ t2 are the 
overall phases of grains 1,2, and o\. 2 are the chiralities 
of the p ± ip order parameter in each grain. 

As we show below, the constants U and J strongly de- 
pend on the behavior of which is in general very dif- 
ficult to obtain from first principles. For a flat interface, 
the component of momentum parallel to the junction is 



conserved, i.e. k 



9\\ 



If the momenta involved are 



close to the Fermi momentum (and assuming for simplic- 
ity a spherically symmetric Fermi surface), then one has 
(approximately) that k 2 + k\ « k 2 F « q 2 + q\\ hence, 
k± « q± or k± « — qx, corresponding to forward and 
backward scattering in the normal direction to the bar- 
rier, respectively. 

There should be a strong suppression of tunneling 
when the vectors k, q are not normal to the interface. 
The reason is that the smaller the perpendicular compo- 
nent, the more exponentially suppressed is the tunneling 
amplitude (for example, consider a WKB approximation: 
the smaller k±, q±, the deeper under the barrier). If S<p is 
a small angle that measures deviations from normal inci- 



dence and 



2" . 



one can show that the main con- 



1 '° "3 

tribution to the Josephson tunneling Hamiltonian comes 
from choosing any of the following four combinations: 



2tt 



2tt 

or <j)f. = —a + Sp + ir (A5) 
o 



7T — 2S(f , 



(A6) 



where the last choice corresponds to forward or backward 
scattering respectively. 

The Josephson coupling can be written in terms of 
these choices as: 

A 2 r 

\t F (5if)\ 2 COS(0 M - 9 2 , a + (oj. - cr 2 ) 6<p) 
-\t B (5ip)\ 2 cos(e ha - 9 2 , a + (<7i + cr 2 ) Sep) , (A7) 

where tp<(5ip), tg(6ip) are forward and backward small 
angle scattering amplitudes (also, recall the definition 
0i M = &i + ^a&i from section Hill. 

Expanding around small Sip before carrying out the 
angular integral, one obtains: 



Ej ~-[U + J trio-a] oos(0i |O - 6 2 ,a) 



(A8) 



Notice that the a-th unit vector e» a is normal to the where 



r = ^- I dS v [\t F (S v )\ 2 -\t B (S v )\ 2 ] (l-Sp) 2 (A9) 



and 



/ = =^ I dSp [MMI 2 + \t B (Sp)\ 2 ] S V 2 . (A10) 



As we discussed above, the barrier is more transpar- 
ent for close to normal incidence, and can be engineered 
so that Sip must remain small, and thus the ratio J/U as 
obtained above can be made controllably small. The pre- 
cise condition for having J <g; U depends on the details 
of tp t B(S(p). As a simple example, for tunneling through 
a square barrier in ordinary quantum mechanics, the ra- 
tio J/U will depend on the height of the barrier V and 
on /cpa, where a is the length of the barrier. The larger 
kpa, the smaller J/U. This model may not capture in 
full detail the underlying physics of the Josephson cou- 
pling problem 2 ^; nevertheless, simple as it is, it shows 
how the structure of the barrier can be used to tune the 
ratio J/U. 

If J <C U, then in the temperature regime J <C T <C U 
the system is effectively constrained to the three-color 
manifold of states: G\. a — 62. a = (mod27r). In this 
case, the effective Hamiltonian for the coupling between 
triangles 1,2 is simply: 



H1.2 = -Jcr 1 a 2 , 
with J > (ferromagnetic coupling). 



(All) 
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